Importing Cell Profiler Images to Python¶
PI Sonja Aits, Lund University¶
Nikolay Oskolkov, SciLifeLab, NBIS Long Term Support, nikolay.oskolkov@scilifelab.se¶
Abstract
Here we will show how to convert Cell Profiler (Cellomics) format C01 to 8-bit tiff images and import those images into Python for further use with scikit-image and Convolutional Neural Networks (CNN). We will also demonstrate basic image processing steps with PIL / Pillow scikit-image Python modules.
Table of Contents:¶
Convert C01 to TIFF ¶
Images which were received from Sonja Aits were in C01 format which is an intrinsic Cell Profiler (Cellomics) format.
Unfortunately, Cell Profiler is primarily a Windows program which is almost uninstallable under Linux. Following Cell Profiler instrictions here https://github.com/CellProfiler/CellProfiler/wiki/Source-Installation-(Ubuntu-14.04-LTS) to install it from source on Ubuntu 14.04 LTS which is my version, I ran into multiple problems with ilastik module https://www.ilastik.org/documentation/basics/installation.html which does not seem to be possible to install under Linux.
from IPython.display import Image
Image('/home/nikolay/WABI/S_Aits/Presentation/CellProfiler2Python/CellProfiler.png')
Therefore I decided to use Bio-Formats library for processing biological microscopy images. More specifically, I used the command line tool bftools which can be downloaded here https://www.openmicroscopy.org/bio-formats/downloads/, bftools includes a very handy tool bfconvert which converts between over 100 formats supported by Bio-Formats, including C01 Cell Profiler format. To run bfconvert let us first list the C01 files and build names for corresponding tiff-files:
import glob
c01_files = glob.glob('/home/nikolay/WABI/S_Aits/co1_images/*.C01')
c01_files
tiff_files = [str(i.split('.C01')[0]) + '.tiff' for i in c01_files]
tiff_files = [i.replace('co1_images','tiff_images') for i in tiff_files]
tiff_files
To run Unix commands such as bfconvert in Python environment we will use subprocess module:
import subprocess
for i,j in zip(c01_files, tiff_files):
print('Converting {} to tiff format'.format(i))
subprocess.run(['/home/nikolay/WABI/S_Aits/bftools/bfconvert', '-overwrite', '-nogroup', i, j])
Now we have converted C01 to tiff-files, however the problem with the tiff-files is that they are 16-bit files which can not be open properly with many graph editors and viewers. We need to use the standard Linux convert command to generate 8-bit tiff-files:
%%bash
for i in $(ls /home/nikolay/WABI/S_Aits/tiff_images/*.tiff)
do
echo Conevrt file $i to 8-bit tiff format
convert $i -auto-level -depth 8 -define quantum:format=unsigned -type grayscale ${i::-5}_8bit.tiff
rm $i
mv ${i::-5}_8bit.tiff $i
done
Hurray! Mission accomplished, we have converted C01-files to tiff format which can be opened by standard Windows / Mac / Linux graph viewers and editors. Now it is time to import those tiff-images into Python using special libraries for image editing.
Reading Images with ImageJ ¶
Cell Profiler is hard to install on Linux but it is possible to install another Bio-Formats module which is called ImageJ. This way one can at least display the C01 images if not convert them to any other useful format. To display the images one needs to perform the following simple operations:
- cd /home/nikolay/Downloads/ImageJ
- ./ImageJ
- Plugins -> Bio-Formats -> Bio-Formats Importer
- Plugins -> Filters -> Color Transformer
from IPython.display import Image
Image('/home/nikolay/WABI/S_Aits/Presentation/CellProfiler2Python/ImageJ.png', width = 2000)
Import Images to Python ¶
Let us now try to import the tiff-images into Python for further image processing with Pillow and scikit-images Python modules. There are two common ways to import images into Python. First, we can use use scikit-image which immediately converts the image into numpy array:
from skimage import io
img = io.imread('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')[0][0][0]
img.shape
img
import matplotlib.pyplot as plt
from skimage import io
fig = plt.figure(figsize=(20, 15))
for i,j in zip(tiff_files,range(len(tiff_files))):
img = io.imread(i)[0][0][0]
plt.subplot(231 + j)
plt.imshow(img[:,:], cmap=plt.cm.gray)
plt.title(i.split('/')[6].split('.')[0], fontsize = 20)
plt.tight_layout()
plt.show()
Superimposing Color Channels ¶
Second, we can use PIL / Pillow module which is also very handy, but here we need to additionally convert the image into numpy array for future usage with Convolutional Neural Networks (CNNs). A great thing about PIL / Pillow module is that one can easily superimpose different color channels onto each other and get a colored image out of 3 black-and-white images. Let us demonstrate how to do it:
import matplotlib.pyplot as plt
from PIL import Image
im1 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')
im2 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d1.tiff')
im3 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d2.tiff')
fig = plt.figure(figsize=(20, 15))
plt.figure(figsize = (12, 8))
plt.subplot(131)
plt.imshow(im1)
plt.subplot(132)
plt.imshow(im2)
plt.subplot(133)
plt.imshow(im3)
plt.tight_layout()
plt.show()
Again, we can easily convert images into numpy arrays:
import numpy as np
imarray1 = np.array(im1)
imarray2 = np.array(im2)
imarray3 = np.array(im3)
imarray1.shape
imarray1
And finally, let us merge black-and-white images from R, G and B color channels and display the result which will be a colored Image.
newImage = Image.merge("RGB", [im1,im2,im3])
np.array(newImage).shape
plt.figure(figsize = (12, 8))
plt.imshow(newImage)
plt.show()
Image Processing ¶
Image processing is the process of creating a new image from an existing image, typically simplifying or enhancing the content in some way. It is a type of digital signal processing and is not concerned with understanding the content of an image. Examples of image processing include:
- Normalizing photometric properties of the image, such as brightness or color.
- Cropping the bounds of the image, such as centering an object in a photograph.
- Removing digital noise from an image, such as digital artifacts from low light levels.
We will perform image processing using Pillow Python module. Pillow is a PIL library that supports Python 3 and is the preferred modern library for image manipulation in Python. It is even required for simple image loading and saving in other Python scientific libraries such as SciPy and Matplotlib.
from PIL import Image
test_image = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')
print(test_image.format)
print(test_image.mode)
print(test_image.size)
import matplotlib.pyplot as plt
plt.imshow(test_image)
plt.show()
Here mode 'L' means that this is an RGB image which was converted to grayscale. Another way to do image processing is in Matplotlib as well but it will run Pillow unders the hood.
# load and display an image with Matplotlib
from matplotlib import image
from matplotlib import pyplot
# load image as pixel array
data = image.imread('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')
# summarize shape of the pixel array
print(data.dtype)
print(data.shape)
# display the array of pixels as an image
pyplot.imshow(data)
pyplot.show()
We can access the pixel data from a Pillow Image. Perhaps the simplest way is to construct a NumPy array and pass in the Image object. The process can be reversed converting a given array of pixel data into a Pillow Image object using the Image.fromarray() function. This can be useful if image data is manipulated as a NumPy array and you then want to save it later as a PNG or JPEG file.
from numpy import asarray
data = asarray(test_image)
# summarize shape
print(data.shape)
# create Pillow image
image2 = Image.fromarray(data)
# summarize image details
print(image2.format)
print(image2.mode)
print(image2.size)
Both approaches are effective for loading image data into NumPy arrays, although the Matplotlib imread() function uses fewer lines of code than loading and converting a Pillow Image object and may be preferred. For example, you could easily load all images in a directory as a list as follows:
# load all images in a directory
from os import listdir
from matplotlib import image
import matplotlib.pyplot as plt
# load all images in a directory
loaded_images = list()
fig = plt.figure(figsize=(20, 15))
tiff_files = listdir('/home/nikolay/WABI/S_Aits/tiff_images/')
for filename,j in zip(tiff_files,range(len(tiff_files))):
# load image
img_data = image.imread('/home/nikolay/WABI/S_Aits/tiff_images/' + filename)
plt.subplot(231 + j)
plt.imshow(img_data)
# store loaded image
loaded_images.append(img_data)
print('> loaded %s %s' % (filename, img_data.shape))
plt.show()
An image object can be saved by calling the save() function. This can be useful if you want to save an image in a different format, in which case the ‘format‘ argument can be specified, such as PNG, GIF, or PEG. Saving images is useful if you perform some data preparation on the image before modeling. One example is converting color images (RGB channels) to grayscale (1 channel). There are a number of ways to convert an image to grayscale, but Pillow provides the convert() function and the mode ‘L‘ will convert an image to grayscale.
print(test_image.mode)
gs_image = test_image.convert(mode='L')
print(gs_image.mode)
np.array(gs_image).shape
plt.imshow(test_image)
plt.show()
plt.imshow(gs_image)
plt.show()
In our case, the images were already in the grayscale mode, so conversion did not change anything, i.e. mode was 'L' before conversion and it remained 'L' after conversion.
Our images are of a very high resolution which might be a problem for the CNN as those images demand a lot of RAM even if one uses minibatches. Here we will show how to downsize all the images and save them with a lower resolution. For example, the test photograph we have been working with has the width and height of (1104, 1104). We can resize it to (100, 100), in which case the largest dimension, e.g. the width, will be reduced to 100, and the height will be scaled in order to retain the aspect ratio of the image. In our case, we have equal width and height so the final dimension will be (100, 100).
import os
import numpy as np
from PIL import Image
tiff_files = os.listdir('/home/nikolay/WABI/S_Aits/tiff_images/')
tiff_downsized_files = os.mkdir('/home/nikolay/WABI/S_Aits/tiff_downsized_images/')
for filename in tiff_files:
img_data = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/' + filename)
print('> loaded %s %s' % (filename, np.array(img_data).shape))
img_data.thumbnail((100,100))
print('> downsized image: {}'.format(np.array(img_data).shape))
print('> saving downsized image {}'.format(filename))
img_data.save('/home/nikolay/WABI/S_Aits/tiff_downsized_images/' + filename)
Simple image manipulation can be used to create new versions of images that, in turn, can provide a richer training dataset when modeling. Generally, this is referred to as data augmentation and may involve creating flipped, rotated, cropped, or other modified versions of the original images with the hope that the algorithm will learn to extract the same features from the image data regardless of where they might appear.
# horizontal flip
hoz_flip = test_image.transpose(Image.FLIP_LEFT_RIGHT)
# vertical flip
ver_flip = test_image.transpose(Image.FLIP_TOP_BOTTOM)
# plot all three images using matplotlib
fig = plt.figure(figsize=(20, 15))
pyplot.subplot(131)
pyplot.title('Original Image', fontsize=20)
pyplot.imshow(test_image)
pyplot.subplot(132)
pyplot.imshow(hoz_flip)
pyplot.title('Horizontal Flip', fontsize=20)
pyplot.subplot(133)
pyplot.imshow(ver_flip)
pyplot.title('Vertical Flip', fontsize=20)
pyplot.show()
Let us now create a simple augmented data set byt flipping each of the 6 downsized images horizontally and vertically:
import os
import numpy as np
from PIL import Image
tiff_files = os.listdir('/home/nikolay/WABI/S_Aits/tiff_downsized_images/')
tiff_augmented_files = os.mkdir('/home/nikolay/WABI/S_Aits/tiff_augmented_images/')
for filename in tiff_files:
img_data = Image.open('/home/nikolay/WABI/S_Aits/tiff_downsized_images/' + filename)
print('> loaded %s %s' % (filename, np.array(img_data).shape))
hoz_flip = img_data.transpose(Image.FLIP_LEFT_RIGHT)
ver_flip = img_data.transpose(Image.FLIP_TOP_BOTTOM)
print('> saving augmented images for {}'.format(filename))
img_data.save('/home/nikolay/WABI/S_Aits/tiff_augmented_images/' + filename)
hoz_flip.save('/home/nikolay/WABI/S_Aits/tiff_augmented_images/aug_hoz_' + filename)
ver_flip.save('/home/nikolay/WABI/S_Aits/tiff_augmented_images/aug_ver_' + filename)
An image can be rotated using the rotate() function and passing in the angle for the rotation. The function offers additional control such as whether or not to expand the dimensions of the image to fit the rotated pixel values (default is to clip to the same size), where to center the rotation the image (default is the center), and the fill color for pixels outside of the image (default is black).
# create rotated versions of an image
from PIL import Image
from matplotlib import pyplot
fig = plt.figure(figsize=(20, 15))
pyplot.subplot(131)
pyplot.title('Original Image', fontsize=20)
pyplot.imshow(test_image)
pyplot.subplot(132)
pyplot.imshow(test_image.rotate(45))
pyplot.title('Rotate 45 Degrees', fontsize=20)
pyplot.subplot(133)
pyplot.imshow(test_image.rotate(90))
pyplot.title('Rotate 90 Degrees', fontsize=20)
pyplot.show()
Let us now make a funny thing. We will create an augmented data set by rotating each image from 0 to 360 degrees by 2 degree step. Therefore each image will be replicated by 180 augmented images. Since we have essentially two photographs, 3 channels each, the first photograph will be represented by 180x3 = 540 images, and the second photograph will be represented by other 540 images. The total data set thus will be 1080 images, where half comes from first photograph and the other half comes from the second photograph. Later we can mark augmented images originated from the first photograph as "cases" and the images from the second photograph as "controls" and run a simple case-control image classification to see if the CNN is capable of discriminating between the two photographs. But firss tlet us create the 1080 imges augmented data set:
import os
import numpy as np
from PIL import Image
tiff_files = os.listdir('/home/nikolay/WABI/S_Aits/tiff_downsized_images/')
tiff_augmented_files = os.mkdir('/home/nikolay/WABI/S_Aits/tiff_augmented_images/')
for filename in tiff_files:
print('Augmenting image {}'.format(filename))
for i in range(0,360,2):
img_data = Image.open('/home/nikolay/WABI/S_Aits/tiff_downsized_images/' + filename)
rot_image = img_data.rotate(i)
rot_image.save('/home/nikolay/WABI/S_Aits/tiff_augmented_images/aug_' + str(i) + '_' + filename)
%%bash
cd /home/nikolay/WABI/S_Aits/tiff_augmented_images/
ls *.tiff | wc -l
We can see that we got 1080 augmented images.
%%bash
cd /home/nikolay/WABI/S_Aits/tiff_augmented_images/
find . -name '*f00*' | wc -l
find . -name '*f01*' | wc -l
Out of them 540 have the 'f00' suffix meaning they belong to the first photograph, and 540 have 'f01' suffix implying they belong to the second photograph. This will be our case-control data set which we will split into trainig and test data sets.
What is CNN and how it works? ¶
Convolutional Neural Network (CNN) is a type of an Artificial Neural Network (ANN) which makes sense to use when dat have a spatial structure. This can be a matrix where columns and rows are ordered and their positions can not be changed, i.e. resampling (permutation, bootstrapping) does not make sense. An image is an example of such an ordered matrix since each pixel has a well-defined neighbour along x- and y-axes. In principle, a text, speach or any Time Series experiment can also be viewed as an ordered matrix but it is a less common practice to use CNN for those types of data. There are three types of layers in a Convolutional Neural Network:
- Convolutional Layers.
- Pooling Layers.
- Fully-Connected Layers.
Fully-Connected layer is a standard Dense layer like in Multilayer Perceptron (MLP), in contrast the other two types of layers are CNN-specific. Both Convolutional and Pooling layers have neurons which are called filters, these filters make some transformation of the input data. The filters can "see" only a limited part of the image, e.g. 5x5 pixel region which is called receptive field. One can think about a receptive field as an input for a neuron. The receptive field window can be "dragged" along horizontal or vertical axis of the image with a stride step and produce multiple receptive fileds. After each filter / neuron was applied to the receptive fields, the output of each neuron will be an ordered matrix which is called a feature map. This terminology is applicable for both Convolutional and Pooling layers, the difference between them is in how texactly the filters / neurons do transformation of the data.
from IPython.display import Image
Image('/home/nikolay/WABI/S_Aits/Presentation/ObjectDetection/CNN.png', width = 2000)
Let us go through those layers one-by-one and explain how each of them works, on the way we will introduce common terminology used in CNNs.
Convolutional Layers¶
Suppose the input dimensions of the image are (1, 100, 100), where the picture has 100x100 pixels resolution and one color channel with values varying from 0 to 255, i.e. it is a grayscale image. Let us apply 100 filters to the image with a 3x3 receptive filed and stride 1, meaning we shift the receptive field window 1 horizontal or vertical position. Each filter applies a linear combintation to the values of the receptive field, like in the figure above. The weights of the filters are learned by minimizing the cost function via back-propogation like in normal ANNs. The output of each neuron of the Convolutional Layer is a feature map with "treansformed" color values for each pixel of the original image, so in total we will have 100 feature maps, each of them has slightly reduced dimensions (something like 95x95) compared to the original image.
Pooling Layers¶
Pooling Layers make sense to apply to the feature maps (output) of the Convolutional Layer. A receptive filed of a Pooling Layer typicaly has dimensions 2x2 applied with a stride of 2 and the filter / neuron simply takes the maximum value out of the 4 digits in the 2x2 receptive filed. The number of neurons / filters in the Pooling Layer is always equal to the number of neurons in the Convolutional Layer. As a result, Pooling Layer produces 100 feature maps, each of them has twice reduced dimension, i.e. 50x50, compared to the original image.
Fully-Connected Layers¶
Finally, after applying Covolutional and Pooling layers sequntially a couple of times, the output is flattened, i.e. reshaped into a one-dimensional array and input into an output layer in a traditional ANN fully-connected manner. The output layer generates predictions and has a softmax or sigmoid (for binary classification problems) activation functions applied to each neuron.
CIFAR10 Image Classification with CNN ¶
Here for warming up purposes we will use a Convolutional Neural Network (CNN) built using Keras for image classification task using CIFAR10 image data set. The CIFAR-10 dataset consists of 60,000 photos divided into 10 classes (hence the name CIFAR-10). Classes include common objects such as airplanes, automobiles, birds, cats and so on. The dataset is split in a standard way, where 50,000 images are used for training a model and the remaining 10,000 for evaluating its performance. The photos are in color with red, green and blue components, but are small measuring 32 by 32 pixel squares.
# Plot ad hoc CIFAR10 instances
from keras.datasets import cifar10
from matplotlib import pyplot
from scipy.misc import toimage
# load data
(X_train, y_train), (X_test, y_test) = cifar10.load_data()
# create a grid of 3x3 images
for i in range(0, 9):
pyplot.subplot(330 + 1 + i)
pyplot.imshow(toimage(X_train[i]))
# show the plot
pyplot.show()
# Simple CNN model for CIFAR-10
import numpy
from keras.datasets import cifar10
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Dropout
from keras.layers import Flatten
from keras.constraints import maxnorm
from keras.optimizers import SGD
from keras.layers.convolutional import Conv2D
from keras.layers.convolutional import MaxPooling2D
from keras.utils import np_utils
from keras import backend as K
K.set_image_dim_ordering('th')
# fix random seed for reproducibility
seed = 7
numpy.random.seed(seed)
# load data
(X_train, y_train), (X_test, y_test) = cifar10.load_data()
# normalize inputs from 0-255 to 0.0-1.0
X_train = X_train.astype('float32')
X_test = X_test.astype('float32')
X_train = X_train / 255.0
X_test = X_test / 255.0
# one hot encode outputs
y_train = np_utils.to_categorical(y_train)
y_test = np_utils.to_categorical(y_test)
num_classes = y_test.shape[1]
# Create the model
model = Sequential()
model.add(Conv2D(32, (3, 3), input_shape=(3, 32, 32), padding='same', activation='relu',
kernel_constraint=maxnorm(3)))
model.add(Dropout(0.2))
model.add(Conv2D(32, (3, 3), activation='relu', padding='same', kernel_constraint=maxnorm(3)))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Flatten())
model.add(Dense(512, activation='relu', kernel_constraint=maxnorm(3)))
model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='softmax'))
# Compile model
epochs = 25
lrate = 0.01
decay = lrate/epochs
sgd = SGD(lr=lrate, momentum=0.9, decay=decay, nesterov=False)
model.compile(loss='categorical_crossentropy', optimizer=sgd, metrics=['accuracy'])
print(model.summary())
# Fit the model
model.fit(X_train, y_train, validation_data=(X_test, y_test), epochs=epochs, batch_size=32)
# Final evaluation of the model
scores = model.evaluate(X_test, y_test, verbose=0)
print("Accuracy: %.2f%%" % (scores[1]*100))
X_train.shape
y_train.shape
Cell Image Classification with CNN ¶
Finally here we will implement a vanilla CNN for classification of biological cell images. First, we need read the images using the Pillow Python module and convert them to numpy arrays. Important to note that we will consider all the 6 images as independent even though it is in fact 2 independent images that have 3 RGB channels each:
from PIL import Image
im1 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')
im2 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d1.tiff')
im3 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d2.tiff')
im4 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d0.tiff')
im5 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d1.tiff')
im6 = Image.open('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f01d2.tiff')
#im1 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f00d0.png')
#im2 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f00d1.png')
#im3 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f00d2.png')
#im4 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f01d0.png')
#im5 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f01d1.png')
#im6 = Image.open('/home/nikolay/WABI/S_Aits/png_images/MFGTMPcx7_170802000001_A01f01d2.png')
import numpy as np
imarray1 = np.array(im1)
imarray2 = np.array(im2)
imarray3 = np.array(im3)
imarray4 = np.array(im4)
imarray5 = np.array(im5)
imarray6 = np.array(im6)
Looking at the dimensions of the arrays we conclude that each of them is a 4-dimensional array. First dimension corresponds to the number of samples, which is 1 for each image; second corresponds to the channel, which is 1 since we treat R, G and B channels as independent; finally, the last two dimensions correspond to the numbers of pixels in each 2D image.
imarray1_f = np.expand_dims(np.expand_dims(imarray1, axis=0), axis=0)
imarray2_f = np.expand_dims(np.expand_dims(imarray2, axis=0), axis=0)
imarray3_f = np.expand_dims(np.expand_dims(imarray3, axis=0), axis=0)
imarray4_f = np.expand_dims(np.expand_dims(imarray4, axis=0), axis=0)
imarray5_f = np.expand_dims(np.expand_dims(imarray5, axis=0), axis=0)
imarray6_f = np.expand_dims(np.expand_dims(imarray6, axis=0), axis=0)
imarray6_f.shape
Here we concatenate 4 images and build a training data set, while the rest two concatenated images will be for testing.
X_image_train = np.concatenate((imarray1_f, imarray2_f, imarray3_f, imarray4_f), axis = 0)
X_image_test = np.concatenate((imarray5_f, imarray6_f), axis = 0)
X_image_train.shape
We assign random binary classes, i.e. 0 and 1, to the training and test data sets.
y_image_train = np.array([0, 0, 1, 1])
y_image_test = np.array([0, 1])
Finally, we will construct the model and normalize the data:
import numpy
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Dropout
from keras.layers import Flatten
from keras.constraints import maxnorm
from keras.optimizers import SGD
from keras.layers.convolutional import Conv2D
from keras.layers.convolutional import MaxPooling2D
from keras.utils import np_utils
from keras import backend as K
K.set_image_dim_ordering('th')
# fix random seed for reproducibility
seed = 7
numpy.random.seed(seed)
# normalize inputs from 0-255 to 0.0-1.0
X_image_train = X_image_train.astype('float32')
X_image_test = X_image_test.astype('float32')
X_image_train = X_image_train / 255.0
X_image_test = X_image_test / 255.0
# one hot encode outputs
y_image_train = np_utils.to_categorical(y_image_train)
y_image_test = np_utils.to_categorical(y_image_test)
num_classes = y_image_test.shape[1]
# Create the model
model = Sequential()
model.add(Conv2D(100, (3, 3), input_shape=(1, 1104, 1104), padding='same', activation='relu',
kernel_constraint=maxnorm(3)))
model.add(Dropout(0.2))
model.add(Conv2D(100, (3, 3), activation='relu', padding='same', kernel_constraint=maxnorm(3)))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Flatten())
model.add(Dense(10, activation='relu', kernel_constraint=maxnorm(3)))
model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='sigmoid'))
# Compile model
epochs = 5
lrate = 0.01
decay = lrate/epochs
sgd = SGD(lr=lrate, momentum=0.9, decay=decay, nesterov=False)
model.compile(loss='binary_crossentropy', optimizer=sgd, metrics=['accuracy'])
print(model.summary())
and start training the model:
# Fit the model
model.fit(X_image_train, y_image_train, validation_data=(X_image_test, y_image_test),
epochs=epochs, batch_size=1)
# Final evaluation of the model
scores = model.evaluate(X_image_test, y_image_test, verbose=0)
print("Accuracy: %.2f%%" % (scores[1]*100))
Very good news here is that the CNN undertands the input and does not throw any errors. The output of course does not make sense since we generated our data sets in a strange way, i.e. assigning random classes to random images. This was done on purpose, here now most important is to check that the input can be understood by the model, i.e. there are no syntax errors, so the mission is accomplished.
Case-Control Cell Image Classification with CNN ¶
Here we will implement a simple Case-Control classification with CNN using the augmented data set, where "cases" will be images from the first photograph and "controls" will be augmented images from the second photograph. we will start with reading images, converting them to 4-dimensional numpy arrays and building X_image and y_image data set, where the former contains the numpy arrays of pixels, while the latter contains labels.
import os
import numpy as np
from PIL import Image
tiff_augmented_files = os.listdir('/home/nikolay/WABI/S_Aits/tiff_augmented_images/')
X_image = np.zeros(shape = (0,1,100,100))
y_image = list()
for filename in tiff_augmented_files:
#print('Working with file {}'.format(filename))
img = Image.open(('/home/nikolay/WABI/S_Aits/tiff_augmented_images/' + filename))
imgarray = np.array(img)
imgarray = np.expand_dims(np.expand_dims(imgarray, axis=0), axis=0)
#print(imgarray.shape)
X_image = np.concatenate((X_image, imgarray), axis = 0)
if('f00' in filename):
y_image.append(0)
else:
y_image.append(1)
print(X_image.shape)
y_image = np.array(y_image)
print(y_image.shape)
print(y_image[0:10])
Let us now split the data set into train and test sub-sets using the scikit-learn handy function:
from sklearn.model_selection import train_test_split
X_image_train,X_image_test, y_image_train,y_image_test = train_test_split(X_image, y_image, test_size = 0.2)
print("The dimensions of X_image_train and y_image_train are",X_image_train.shape,"\b and",
y_image_train.shape,"\b, respectively")
print("The dimensions of X_image_test and y_image_test are",
X_image_test.shape,"\b and",y_image_test.shape,"\b, respectively")
Now let us build a simple CNN model where we will use data normalization (by max value) and one-hot encoding of the labels:
# Simple CNN model for CIFAR-10
import numpy
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Dropout
from keras.layers import Flatten
from keras.constraints import maxnorm
from keras.optimizers import SGD
from keras.layers.convolutional import Conv2D
from keras.layers.convolutional import MaxPooling2D
from keras.utils import np_utils
from keras import backend as K
K.set_image_dim_ordering('th')
# fix random seed for reproducibility
seed = 7
numpy.random.seed(seed)
# normalize inputs from 0-255 to 0.0-1.0
X_image_train = X_image_train.astype('float32')
X_image_test = X_image_test.astype('float32')
X_image_train = X_image_train / 255.0
X_image_test = X_image_test / 255.0
# one hot encode outputs
y_image_train = np_utils.to_categorical(y_image_train)
y_image_test = np_utils.to_categorical(y_image_test)
num_classes = y_image_test.shape[1]
# Create the model
model = Sequential()
model.add(Conv2D(100, kernel_size=(3, 3), input_shape=(1, 100, 100), padding='same', activation='relu',
kernel_constraint=maxnorm(3)))
model.add(Dropout(0.2))
model.add(Conv2D(100, kernel_size=(3, 3), activation='relu', padding='same', kernel_constraint=maxnorm(3)))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Flatten())
model.add(Dense(512, activation = 'relu', kernel_constraint = maxnorm(3)))
model.add(Dropout(0.5))
model.add(Dense(num_classes, activation = 'sigmoid'))
# Compile model
epochs = 5
lrate = 0.01
decay = lrate / epochs
sgd = SGD(lr = lrate, momentum = 0.9, decay = decay, nesterov = False)
model.compile(loss = 'binary_crossentropy', optimizer = sgd, metrics = ['accuracy'])
print(model.summary())
import pydot_ng as pydot
import graphviz
from IPython.display import SVG
from keras.utils.vis_utils import model_to_dot
SVG(model_to_dot(model).create(prog='dot', format='svg'))
# Fit the model
estimator = model.fit(X_image_train, y_image_train, validation_split = 0.2, shuffle = True,
epochs = epochs, batch_size = 5)
print("Training Accuracy: ",estimator.history['acc'][-1])
print("Validation Accuracy: ",estimator.history['val_acc'][-1])
# summarize history for accuracy
plt.plot(estimator.history['acc'])
plt.plot(estimator.history['val_acc'])
plt.title('Model Accuracy', fontsize = 20)
plt.ylabel('Accuracy', fontsize = 20)
plt.xlabel('Epoch', fontsize = 20)
plt.legend(['Train', 'Validate'], loc='upper left', fontsize = 20)
plt.show()
# summarize history for loss
plt.plot(estimator.history['loss'])
plt.plot(estimator.history['val_loss'])
plt.title('Model Loss', fontsize = 20)
plt.ylabel('Loss', fontsize = 20)
plt.xlabel('Epoch', fontsize = 20)
plt.legend(['Train', 'Validate'], loc='upper left', fontsize = 20)
plt.show()
We observe a very fast (after only 5 epochs) decrease of the loss function and the accuracy quickly reaches 100% on the validation data set.
# Final evaluation of the model
scores = model.evaluate(X_image_test, y_image_test, verbose=0)
print("Accuracy: %.2f%%" % (scores[1]*100))
ann_pred_probs = model.predict(X_image_test)
ann_pred_probs[0:5]
y_image_test[0:5]
Thus we reach a very high accuracy of discriminating between first and second photographs, which is to be expected since those two photograps visually look very different, still good to see the demonstration of Computer Vision. We can see that the predicted labels of the images almost always agree with the true labels.
Object Detection ¶
Images received from Sonja Aits were in C01 format which is an intrinsic Cell Profiler (Cellomics) format. We did image processing and augmentation using Pillow Python module. Now we will implement OpenCV Python module in order to detect cells on the images.
import cv2
import numpy as np
import matplotlib.pyplot as plt
img = cv2.imread('/home/nikolay/WABI/S_Aits/tiff_images/MFGTMPcx7_170802000001_A01f00d0.tiff')
plt.imshow(img)
plt.show()
We are going to do some simple image manipulation: turn the image to grayscale, binarize and dilate using custom kernels.
fig = plt.figure(figsize=(20, 15))
#grayscale
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
plt.subplot(131)
plt.imshow(gray)
plt.title('Grayscale')
#binary
ret,thresh = cv2.threshold(gray,127,255,0)
plt.subplot(132)
plt.imshow(thresh)
plt.title('Binary')
#dilation
kernel = np.ones((1,1), np.uint8)
img_dilation = cv2.dilate(thresh, kernel, iterations=1)
plt.subplot(133)
plt.imshow(img_dilation)
plt.title('Dilated')
plt.show()
Next, we will find contours using the OpenCV function findCoutours:
# find contours
ctrs, hier = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
# sort contours
sorted_ctrs = sorted(ctrs, key=lambda ctr: cv2.boundingRect(ctr)[0])
cv2.drawContours(img, ctrs, -1, (0,255,0), 3)
for i, ctr in enumerate(sorted_ctrs):
# Get bounding box
x, y, w, h = cv2.boundingRect(ctr)
# Getting ROI
roi = img[y:y + h, x:x + w]
# show ROI
# cv2.imshow('segment no:'+str(i),roi)
cv2.rectangle(img, (x, y), (x + w, y + h), (0, 255, 0), 2)
if w > 15 and h > 15:
cv2.imwrite('/home/nikolay/WABI/S_Aits/bounding_boxes/i.png'.format(i), roi)
fig = plt.figure(figsize=(20, 15))
plt.imshow(img)
plt.show()
Now we have detected the cell on the images, put bounding boxes around them, croped the cell amd written them to separate imge files. Next step is to annotate the bounding boxes and input them to Faster RCNN.